The relative benefits of environmental enrichment on learning and memory are greater when stressed: a meta-analysis of interactions in rodents
Supplementary Material S3
Setting-up
Loading packages
# ggthemr needs to be dowlonaded from a github repo
# devtools::install_github('Mikata-Project/ggthemr', force = TRUE)
pacman::p_load(tidyverse, here, metafor, clubSandwich, orchaRd, MuMIn, patchwork,
GoodmanKruskal, networkD3, ggplot2, ggsignif, visdat, ggalluvial, ggthemr, cowplot,
grDevices, png, grid, gridGraphics, pander, formatR)
# needed for model selection using MuMIn within metafor
eval(metafor:::.MuMIn)Loading data and functions
This loads the unprocessed data file and custom functions including
- functions for calculating ‘focal’ effect sizes and variance for lnRR and SMD:
effect_set - functions for calculating ‘pairwise’ comparisons (lnRR) and variance:
effect_set2 - functions for calculating pair-wise contrasts between moderators:
contrast_fun
dat <- read_csv(here("Data", "Data_raw.csv"))
# Load custom function to extract data
source(here("R/Functions.R"))Data organisation
- removing one study (Wang et al. 2020) with negative values as lnRR cannot be calculated with negative values
- rounding down sample sizes that are coded as decimals due to averaging sample sizes across treatments
- calculating effect sizes using custom functions
- ‘flipping’ effect sizes so that all effect sizes that are higher values can be interpreted as individuals performing better in cognitive assays
- assigning human readable terms to moderators initially coded as numbers
- creating a variance-covariance matrix (VCV)
# removing study with negative values as these are unable to be used for lnRR
dat <- droplevels(dat[!dat$First_author == "Wang", ])
# rounding down sample sizes
dat$CC_n <- floor(dat$CC_n)
dat$EC_n <- floor(dat$EC_n)
dat$CS_n <- floor(dat$CS_n)
dat$ES_n <- floor(dat$CS_n)
# 'Focal' effect_size
effect_size <- with(dat, mapply(effect_set, CC_n, CC_mean, CC_SD, EC_n, EC_mean,
EC_SD, CS_n, CS_mean, CS_SD, ES_n, ES_mean, ES_SD, percent = Response_percent,
SIMPLIFY = FALSE))
effect_size <- map_dfr(effect_size, I)
# 'Pairwise' effect size
effect_size2 <- with(dat, mapply(effect_set2, CC_n, CC_mean, CC_SD, EC_n, EC_mean,
EC_SD, CS_n, CS_mean, CS_SD, ES_n, ES_mean, ES_SD, percent = Response_percent,
SIMPLIFY = FALSE))
effect_size2 <- map_dfr(effect_size2, I)
full_info <- which(complete.cases(effect_size) == TRUE)
# adding effect sizes as column
dat <- bind_cols(dat, effect_size, effect_size2)
dat <- dat[full_info, ]
# Flipping 'lower is better' to 'higher is better' effect sizes flipping lnRR
# for values where higher = worse
dat$lnRR_Ea <- ifelse(dat$Response_direction == 2, dat$lnRR_E * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$lnRR_E))
# currently NAswhich causes error
dat$lnRR_Sa <- ifelse(dat$Response_direction == 2, dat$lnRR_S * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$lnRR_S)) # currently NAswhich causes error
dat$lnRR_ESa <- ifelse(dat$Response_direction == 2, dat$lnRR_ES * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$lnRR_ES)) # currently NAswhich causes error
# flipping 'pure effect sizes'
dat$lnRR_E2a <- ifelse(dat$Response_direction == 2, dat$lnRR_E2 * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$lnRR_E2)) # currently NAswhich causes error
dat$lnRR_S2a <- ifelse(dat$Response_direction == 2, dat$lnRR_S2 * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$lnRR_S2)) # currently NAswhich causes error
dat$lnRR_ES2a <- ifelse(dat$Response_direction == 2, dat$lnRR_ES2 * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$lnRR_ES2)) # currently NAswhich causes error
dat$lnRR_E3a <- ifelse(dat$Response_direction == 2, dat$lnRR_E3 * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$lnRR_E3)) # currently NAswhich causes error
dat$lnRR_S3a <- ifelse(dat$Response_direction == 2, dat$lnRR_S3 * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$lnRR_S3)) # currently NAswhich causes error
# flipping SMD
dat$SMD_Ea <- ifelse(dat$Response_direction == 2, dat$SMD_E * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$SMD_E)) # currently NAswhich causes error
dat$SMD_Sa <- ifelse(dat$Response_direction == 2, dat$SMD_S * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$SMD_S)) # currently NAswhich causes error
dat$SMD_ESa <- ifelse(dat$Response_direction == 2, dat$SMD_ES * -1, ifelse(is.na(dat$Response_direction) ==
TRUE, NA, dat$SMD_ES))
# assigning human readable term
dat <- dat %>%
mutate(Type_assay = case_when(Type_assay == 1 ~ "Habituation", Type_assay ==
2 ~ "Conditioning", Type_assay == 3 ~ "Recognition", Type_assay == 4 ~ "Unclear"),
Learning_vs_memory = case_when(Learning_vs_memory == 1 ~ "Learning", Learning_vs_memory ==
2 ~ "Memory", Learning_vs_memory == 3 ~ "Habituation"), Type_reinforcement = case_when(Type_reinforcement ==
1 ~ "Appetitive", Type_reinforcement == 2 ~ "Aversive", Type_reinforcement ==
3 ~ "Not applicable", Type_reinforcement == 4 ~ "Unclear"), Type_stress_exposure = case_when(Type_stress_exposure ==
1 ~ "Density", Type_stress_exposure == 2 ~ "Scent", Type_stress_exposure ==
3 ~ "Shock", Type_stress_exposure == 4 ~ "Exertion", Type_stress_exposure ==
5 ~ "Restraint", Type_stress_exposure == 6 ~ "Maternal separation", Type_stress_exposure ==
7 ~ "Circadian rhythm", Type_stress_exposure == 8 ~ "Noise", Type_stress_exposure ==
9 ~ "Other", Type_stress_exposure == 10 ~ "Combination", Type_stress_exposure ==
11 ~ "unclear"), Age_stress_exposure = case_when(Age_stress_exposure ==
1 ~ "Prenatal", Age_stress_exposure == 2 ~ "Early postnatal", Age_stress_exposure ==
3 ~ "Adolescent", Age_stress_exposure == 4 ~ "Adult", Age_stress_exposure ==
5 ~ "Unclear"), Stress_duration = case_when(Stress_duration == 1 ~ "Acute",
Stress_duration == 2 ~ "Chronic", Stress_duration == 3 ~ "Intermittent",
Stress_duration == 4 ~ "Unclear"), EE_social = case_when(EE_social ==
1 ~ "Social", EE_social == 2 ~ "Non-social", EE_social == 3 ~ "Unclear"),
EE_exercise = case_when(EE_exercise == 1 ~ "Exercise", EE_exercise == 2 ~
"No exercise"), Age_EE_exposure = case_when(Age_EE_exposure == 1 ~ "Prenatal",
Age_EE_exposure == 2 ~ "Early postnatal", Age_EE_exposure == 3 ~ "Adolescent",
Age_EE_exposure == 4 ~ "Adult", Age_EE_exposure == 5 ~ "Unclear"), Exposure_order = case_when(Exposure_order ==
1 ~ "Stress first", Exposure_order == 2 ~ "Enrichment first", Exposure_order ==
3 ~ "Concurrently", Exposure_order == 4 ~ "Unclear"), Age_assay = case_when(Age_assay ==
1 ~ "Early postnatal", Age_assay == 2 ~ "Adolescent", Age_assay == 3 ~
"Adult", Age_assay == 4 ~ "Unclear"), Sex = case_when(Sex == 1 ~ "Female",
Sex == 2 ~ "Male", Sex == 3 ~ "Mixed", Sex == 4 ~ "Unclear"), Type_EE_exposure = case_when(Type_EE_exposure ==
1 ~ "Nesting material", Type_EE_exposure == 2 ~ "Objects", Type_EE_exposure ==
3 ~ "Cage complexity", Type_EE_exposure == 4 ~ "Wheel/trademill", Type_EE_exposure ==
5 ~ "Combination", Type_EE_exposure == 6 ~ "Other", Type_EE_exposure ==
7 ~ "Unclear"), ROB_blinding = case_when(ROB_blinding == 1 ~ "Yes", ROB_blinding ==
2 ~ "No", ROB_blinding == 3 ~ "Unclear"), ROB_randomisation = case_when(ROB_randomisation ==
1 ~ "Yes", ROB_randomisation == 2 ~ "No", ROB_randomisation == 3 ~ "Unclear"))
# making variance VCVs
VCV_E <- impute_covariance_matrix(vi = dat$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
VCV_S <- impute_covariance_matrix(vi = dat$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ES <- impute_covariance_matrix(vi = dat$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
VCV_Ea <- impute_covariance_matrix(vi = dat$SMDV_E, cluster = dat$Study_ID, r = 0.5)
VCV_Sa <- impute_covariance_matrix(vi = dat$SMDV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ESa <- impute_covariance_matrix(vi = dat$SMDV_ES, cluster = dat$Study_ID, r = 0.5)
# write.csv(dat, file = here('Data', 'Data_processed.csv'), row.names = TRUE)Data exploration
General
Number of effect sizes:
# Number of effect sizes
length(unique(dat$ES_ID))## [1] 92
Number of studies:
length(unique(dat$Study_ID))## [1] 30
Publication year range:
min(dat$Year_published)## [1] 2006
max(dat$Year_published)## [1] 2021
Visual of missing data
# only selecting the columns used in the analysis
dat_missing <- dat %>%
select(c(Study_ID, Year_published, Common_species, Strain, Sex, Age_EE_exposure,
Age_stress_exposure, Age_assay, EE_exercise, EE_social, Type_stress_exposure,
Stress_duration, Exposure_order, Learning_vs_memory, Type_assay, Type_reinforcement,
lnRR_E, lnRRV_E, lnRR_S, lnRRV_S, lnRR_ES, lnRRV_ES, lnRR_E2, lnRRV_E2, lnRR_S2,
lnRRV_S2, lnRR_ES2, lnRRV_ES2, SMD_E, SMDV_E, SMD_S, SMDV_S, SMD_ES, SMDV_ES))
plot_missing <- vis_miss(dat_missing) + theme(plot.title = element_text(hjust = 0.5,
vjust = 3), plot.margin = margin(t = 0.5, r = 3, b = 1, l = 1, unit = "cm")) +
ggtitle("Missing data in the selected predictors") #no missing values
plot_missing# useGoodman and Kruskal’s τ measure of association between categorical
# predictor variables (function from package GoodmanKruskal:
# https://cran.r-project.org/web/packages/GoodmanKruskal/vignettes/GoodmanKruskal.html)
# GKmatrix <- GKtauDataframe(subset(dat, select = c('Sex', 'Type_assay',
# 'Learning_vs_memory', #'Type_reinforcement', 'Type_stress_exposure',
# 'Age_stress_exposure', 'Stress_duration', #'EE_social_HR', 'EE_exercise',
# 'Age_EE_exposure', 'Exposure_order', 'Age_assay'))) plot(GKmatrix)
# simple pairwise contingency tables table(dat$Type_assay,
# dat$Type_reinforcement) table(dat$Age_stress_exposure, dat$Age_EE_exposure)
# table(dat$Type_stress_exposure, dat$Age_stress_exposure)
# table(dat$Type_stress_exposure, dat$Age_assay)
# table(dat$Type_stress_exposure, dat$Stress_duration)Fig. S1 Visual of missing data for each column. No missing data for moderators, random effects, or effect sizes.
Alluvial diagrams
Shows the relationship/nestedness of different data elements (used to produce Fig. 2)
Subjects info: species-strain-sex
Figure corresponds to Fig. 2A in main text
freq_A <- as.data.frame(table(dat$Sex, dat$Common_species, dat$Strain)) %>%
rename(Sex = Var1, Species = Var2, Strain = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_A), axes = 1:3, silent = TRUE)
p1 <- ggplot(data = freq_A, aes(axis1 = Sex, axis2 = Species, axis3 = Strain, y = Freq)) +
geom_alluvium(aes(fill = Sex)) + geom_flow() + geom_stratum(aes(fill = Sex)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + geom_text(stat
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + "stratum",
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + aes(label
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + after_stat(stratum)))
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + #theme_minimal()
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
theme_void() + theme(legend.position = "none", plot.title = element_text(hjust = 0,
vjust = 3), axis.title.x = element_text(), axis.text.x = element_text(face = "bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) + scale_x_discrete(limits = c("Sex",
"Species", "Strain"), expand = c(0.15, 0.05), position = "top") + ggtitle("A study subjects")
p1Fig. S2a Relationship/nestedness between study species, strain of rodent, and sex
Stress info: age-duration-type stress
Figure corresponds to Fig. 2C in main text
freq_C <- as.data.frame(table(dat$Age_stress_exposure, dat$Stress_duration, dat$Type_stress_exposure)) %>%
rename(Age_stress = Var1, Duration_stress = Var2, Type_stress = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_C), axes = 1:3, silent = TRUE)
p3 <- ggplot(data = freq_C, aes(axis1 = Age_stress, axis2 = Duration_stress, axis3 = Type_stress,
y = Freq)) + geom_alluvium(aes(fill = Age_stress)) + geom_flow() + geom_stratum(aes(fill = Age_stress)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + geom_text(stat
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + "stratum",
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + aes(label
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + after_stat(stratum)))
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + #theme_minimal()
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
theme_void() + theme(legend.position = "none", plot.title = element_text(hjust = 0,
vjust = 3), axis.title.x = element_text(), axis.text.x = element_text(face = "bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) + scale_x_discrete(limits = c("Age",
"Duration", "Type"), expand = c(0.1, 0.1), position = "top") + ggtitle("C stress exposure")
p3Fig. S2c Relationship/nestedess between the age at stress exposure, the duration of the stress, and the type of stress
Assay info: LvsM-type-reinforcement
Figure corresponds to Fig. S2D in main text
freq_D <- as.data.frame(table(dat$Learning_vs_memory, dat$Type_assay, dat$Type_reinforcement)) %>%
rename(Learning_Memory = Var1, Type = Var2, Reinforcement = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_D), axes = 1:3, silent = TRUE)
p4 <- ggplot(data = freq_D, aes(axis1 = Learning_Memory, axis2 = Type, axis3 = Reinforcement,
y = Freq)) + geom_alluvium(aes(fill = Learning_Memory)) + geom_flow() + geom_stratum(aes(fill = Learning_Memory)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + geom_text(stat
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + "stratum",
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + aes(label
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + after_stat(stratum)))
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + #theme_minimal()
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
theme_void() + theme(legend.position = "none", plot.title = element_text(hjust = 0,
vjust = 3), axis.title.x = element_text(), axis.text.x = element_text(face = "bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) + scale_x_discrete(limits = c("Learning_Memory",
"Type", "Reinforcement"), expand = c(0.1, 0.1), position = "top") + ggtitle("D cognitive assay")
p4Fig. S2d Relationship/nestedness between if the assay measured a learning or memory trait, the type of assay used, and the type of reinforcement used
Combined plot
Corresponds to Fig. 2 in main text
# p1 + scale_fill_brewer(palette = 'Set3') #Pastel1
(p1 + scale_fill_brewer(palette = "Set3"))/(p2 + scale_fill_brewer(palette = "Set3"))/(p3 +
scale_fill_brewer(palette = "Set3"))/(p4 + scale_fill_brewer(palette = "Set3")) +
plot_layout(ncol = 1, heights = c(1, 1, 1, 1, 1))# ggsave(file = './figs/Alluvial_diagram_v2.pdf', width = 10, height = 12,
# units = 'cm', dpi = 300, scale = 2) #, device = cairo_pdf)Fig. S2 The combined plot of all alluvial diagram used to make Fig. 2 in the main text
Risk of Bias
Number of studies that used randomisation:
dat %>%
group_by(ROB_randomisation) %>%
summarise(n = n_distinct(Study_ID))## # A tibble: 2 x 2
## ROB_randomisation n
## <chr> <int>
## 1 Unclear 15
## 2 Yes 15
15/30 #proportion## [1] 0.5
number of studies that used blinding:
# blinding
dat %>%
group_by(ROB_blinding) %>%
summarise(n = n_distinct(Study_ID))## # A tibble: 2 x 2
## ROB_blinding n
## <chr> <int>
## 1 Unclear 24
## 2 Yes 6
6/30 #proportion## [1] 0.2
Modelling with lnRR
Main effect: environmental enrichment
Meta-analysis
To quantify differences in individual performance in cognitive assays with environmental enrichment, we used the logarithm of response ratio (lnRR) calculated as:
\[ \ln{\text{RR}_\text{EE}} = \ln\left( \frac{ M_\text{ES} + M_\text{EC}} {2} \right) - \ln\left( \frac {M_\text{CS} + M_\text{CC}} {2} \right) = \ln \left( {\frac {M_\text{ES} + M_\text{EC}} {M_\text{CS} + M_\text{CC}}} \right), \]
Variance was calculated as:
\[ \text{var}(\ln{\text{RR}_\text{EE}}) = \left( \frac{1} { M_\text{ES} + M_\text{EC} } \right)^2 \left( \frac{SD_\text{ES}^2}{N_\text{ES}} + \frac{SD_\text{EC}^2}{N_\text{EC}} \right) + \left( \frac {1} { M_\text{CS} + M_\text{CC} } \right)^2 \left(\frac{SD_\text{CS}^2}{N_\text{CS}} + \frac{SD_\text{CC}^2}{n_\text{CC}} \right) \] Variance was converted in to variance-covariance (VCV) matrix (see above in ‘Data organisation’) to control for non-independence in sampling variances from the same studies.
# dat <- read_csv(here('Data','Data_processed.csv'))
mod_E0 <- rma.mv(yi = lnRR_Ea, V = VCV_E, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_E0)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.7836 19.5672 27.5672 37.6106 28.0323
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0082 0.0905 30 no Study_ID
## sigma^2.2 0.0345 0.1858 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 809.2712, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1860 0.0320 5.8116 91 <.0001 0.1224 0.2496 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.310673e-01 1.786605e-01 7.524068e-01 6.097041e-10
Figure corresponds to Fig. 3A in main text
orchard_plot(mod_E0, mod = "Int", xlab = "lnRR", alpha=0.4) + # Orchard plot
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5)+ # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2)+ # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_colour_manual(values = "darkorange")+ # change colours
scale_fill_manual(values="darkorange")+
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Fig. S3 Orchard plot showing meta-analytic mean, 95% confidence interval (thick black line), and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Meta-regression: uni-moderator
To explain some of the unexplained variation in the main effect of environmental enrichment meta-analytic model, we conducted a series of uni-moderator analyses. We calculated marginal R^2 for each moderator, as well as conducted a series of pair-wise contrasts between moderator categories.
Any moderator categories with k < 5 were removed.
Learning vs Memory
Was the measured response interpreted as learning or memory?
mod_E1 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Learning_vs_memory - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_E1)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.3793 18.7586 28.7586 41.2577 29.4729
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0076 0.0873 30 no Study_ID
## sigma^2.2 0.0340 0.1843 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 802.5794, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 18.2861, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.2303 0.0450 5.1227 90 <.0001 0.1410
## Learning_vs_memoryMemory 0.1624 0.0355 4.5713 90 <.0001 0.0919
## ci.ub
## Learning_vs_memoryLearning 0.3197 ***
## Learning_vs_memoryMemory 0.2330 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E1)## R2_marginal R2_coditional
## 0.02572583 1.00000000
Figure corresponds to Fig. 4A in main text
LvsM_E <- orchard_plot(mod_E1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21) + # mean estimate
scale_size_continuous(range = c(1, 7)) + # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_EFig. S4a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_E1 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = Learning_vs_memory,
VCV = VCV_E)
res_table_mod_E1 <- get_estimate(model = mod_E1, contra = contra_mod_E1, moderator = Learning_vs_memory)
res_table_mod_E1## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Learning 0.230 0.141 0.320 0.00000170 -0.185 0.645
## 2 Memory 0.162 0.0919 0.233 0.0000154 -0.249 0.574
## 3 Learning-Memory -0.0679 -0.164 0.0284 0.165 -0.484 0.349
Type of assay
The broad category of the type of assay used to measure learning or memory.
dat1 <- filter(dat, Type_assay %in% c("Recognition", "Habituation", "Conditioning"))
VCV_E1 <- impute_covariance_matrix(vi = dat1$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
mod_E2 <- rma.mv(yi = lnRR_Ea, V = VCV_E1, mod = ~Type_assay - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat1)
summary(mod_E2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.8496 15.6993 27.6993 42.6311 28.7237
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0099 0.0995 30 no Study_ID
## sigma^2.2 0.0321 0.1792 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 661.8903, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.7993, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning 0.2167 0.0351 6.1650 89 <.0001 0.1468 0.2865
## Type_assayHabituation 0.1821 0.1147 1.5886 89 0.1157 -0.0457 0.4100
## Type_assayRecognition 0.0554 0.0640 0.8659 89 0.3889 -0.0717 0.1826
##
## Type_assayConditioning ***
## Type_assayHabituation
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E2)## R2_marginal R2_coditional
## 0.07376134 1.00000000
Figure corresponds to Fig. 4B in main text
Learning_E <- orchard_plot(mod_E2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34", "grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Conditioning", "Recognition")),
annotations = "*") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_EFig. S4b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_E2 <- contrast_fun(data = dat1, response = lnRR_Ea, moderator = Type_assay,
VCV = VCV_E1)
res_table_mod_E2 <- get_estimate(model = mod_E2, contra = contra_mod_E2, moderator = Type_assay)
res_table_mod_E2## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Conditioning 0.217 0.147 0.287 2.02e-8 -0.197 0.630
## 2 Habituation 0.182 -0.0457 0.410 1.16e-1 -0.285 0.649
## 3 Recognition 0.0554 -0.0717 0.183 3.89e-1 -0.371 0.482
## 4 Conditioning-Habitua… -0.0345 -0.262 0.193 7.63e-1 -0.501 0.432
## 5 Conditioning-Recogni… -0.161 -0.296 -0.0265 1.96e-2 -0.590 0.268
## 6 Habituation-Recognit… -0.127 -0.381 0.127 3.24e-1 -0.607 0.353
Type of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
dat2 <- filter(dat, Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"))
VCV_E2 <- impute_covariance_matrix(vi = dat2$lnRRV_E, cluster = dat2$Study_ID, r = 0.5)
mod_E3 <- rma.mv(yi = lnRR_Ea, V = VCV_E2, mod = ~Type_reinforcement - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat2)
summary(mod_E3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.8702 15.7405 27.7405 42.6723 28.7649
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0109 0.1046 30 no Study_ID
## sigma^2.2 0.0319 0.1787 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 764.5984, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.4138, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.2175 0.0726 2.9942 89 0.0036 0.0732
## Type_reinforcementAversive 0.2191 0.0410 5.3456 89 <.0001 0.1377
## Type_reinforcementNot applicable 0.0812 0.0560 1.4504 89 0.1505 -0.0300
## ci.ub
## Type_reinforcementAppetitive 0.3618 **
## Type_reinforcementAversive 0.3006 ***
## Type_reinforcementNot applicable 0.1924
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E3)## R2_marginal R2_coditional
## 0.07720103 1.00000000
Figure corresponds to Fig. 4C in main text
Reinforcement_E <-orchard_plot(mod_E3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34", "grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Aversive", "Not applicable")),
annotations = "*") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_EFig. S4c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_E3 <- contrast_fun(data = dat2, response = lnRR_Ea, moderator = Type_reinforcement,
VCV = VCV_E2)
res_table_mod_E3 <- get_estimate(model = mod_E3, contra = contra_mod_E3, moderator = Type_reinforcement)
res_table_mod_E3## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Appetitive 0.217 0.0732 0.362 3.56e-3 -0.218 0.653
## 2 Aversive 0.219 0.138 0.301 6.89e-7 -0.200 0.638
## 3 Not applicable 0.0812 -0.0300 0.192 1.50e-1 -0.345 0.507
## 4 Appetitive-Aversive 0.00164 -0.163 0.167 9.84e-1 -0.442 0.445
## 5 Appetitive-Not applic… -0.136 -0.315 0.0425 1.33e-1 -0.585 0.312
## 6 Aversive-Not applicab… -0.138 -0.259 -0.0172 2.56e-2 -0.567 0.291
Age at enrichment
The effect of age when individuals were exposed to environmental enrichment.
dat3 <- filter(dat, Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_E3 <- impute_covariance_matrix(vi = dat3$lnRRV_E, cluster = dat3$Study_ID, r = 0.5)
mod_E4 <- rma.mv(yi = lnRR_Ea, V = VCV_E3, mod = ~Age_EE_exposure - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat3)
summary(mod_E4)##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -6.9490 13.8980 23.8980 36.1698 24.6480
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0078 0.0885 29 no Study_ID
## sigma^2.2 0.0327 0.1808 88 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 782.6092, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 18.9060, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Age_EE_exposureAdolescent 0.1799 0.0370 4.8553 86 <.0001 0.1062 0.2535
## Age_EE_exposureAdult 0.2340 0.0620 3.7734 86 0.0003 0.1107 0.3573
##
## Age_EE_exposureAdolescent ***
## Age_EE_exposureAdult ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E4)## R2_marginal R2_coditional
## 0.01127347 1.00000000
Figure corresponds to Fig. 4D in main text
Age_E<- orchard_plot(mod_E4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_EFig. S4d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_E4 <- contrast_fun(data = dat3, response = lnRR_Ea, moderator = Age_EE_exposure,
VCV = VCV_E3)
res_table_mod_E4 <- get_estimate(model = mod_E4, contra = contra_mod_E4, moderator = Age_EE_exposure)
res_table_mod_E4## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adolescent 0.180 0.106 0.254 0.00000533 -0.227 0.587
## 2 Adult 0.234 0.111 0.357 0.000295 -0.185 0.653
## 3 Adolescent-Adult 0.0541 -0.0895 0.198 0.456 -0.371 0.479
Exercise enrichment
Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?
mod_E5 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~EE_exercise - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_E5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.0303 20.0606 30.0606 42.5597 30.7749
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0096 0.0980 30 no Study_ID
## sigma^2.2 0.0345 0.1857 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 807.7169, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 16.1666, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1849 0.0407 4.5464 90 <.0001 0.1041 0.2657
## EE_exerciseNo exercise 0.1900 0.0556 3.4151 90 0.0010 0.0795 0.3005
##
## EE_exerciseExercise ***
## EE_exerciseNo exercise ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E5)## R2_marginal R2_coditional
## 0.0001360993 0.9999999987
Figure corresponds to Fig. 4E in main text
Exercise_E <-orchard_plot(mod_E5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Exercise_EFig. S4e Orchard plot showing the group-wise means of the categorical variable ‘EE_exercise’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_E5 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = EE_exercise,
VCV = VCV_E)
res_table_mod_E5 <- get_estimate(model = mod_E5, contra = contra_mod_E5, moderator = EE_exercise)
res_table_mod_E5## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Exercise 0.185 0.104 0.266 0.0000169 -0.240 0.610
## 2 No exercise 0.190 0.0795 0.301 0.000958 -0.241 0.621
## 3 Exercise-No exercise 0.00511 -0.132 0.142 0.941 -0.434 0.444
Publication bias & sensitivity analysis
Multi-moderator model
We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.
# filter data so that all K < 5 are removed
dat_Efm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
c("Appetitive", "Aversive", "Not applicable"), EE_social %in% c("Social",
"Non-social"), Age_EE_exposure %in% c("Adult", "Adolescent"))
# VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster =
# dat_Efm$Study_ID, r = 0.5)
mod_Efm <- rma.mv(yi = lnRR_Sa, V = VCV_Efm, mod = ~Type_assay - 1 + Learning_vs_memory +
Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Efm)
# summary(mod_Efm) r2_ml(mod_Efm)
res_Efm <- dredge(mod_Efm, trace = 2)
saveRDS(res_Efm, file = here("Rdata", "res_Efm.rds"))
# also saving the full model and data
saveRDS(mod_Efm, file = here("Rdata", "mod_Efm.rds"))
saveRDS(dat_Efm, file = here("Rdata", "dat_Efm.rds"))The akaike weights for the top set of models with AIC < 6
dat_Efm <- readRDS(file = here("Rdata", "dat_Efm.rds"))
mod_Efm <- readRDS(file = here("Rdata", "mod_Efm.rds"))
res_Efm <- readRDS(file = here("Rdata", "res_Efm.rds"))
res_Efm2 <- subset(res_Efm, delta <= 6, recalc.weights = FALSE)
importance(res_Efm2)## Age_EE_exposure Type_assay EE_exercise EE_social
## Sum of weights: 0.63 0.36 0.16 0.14
## N containing models: 11 10 5 5
## Learning_vs_memory Type_reinforcement
## Sum of weights: 0.10 0.07
## N containing models: 4 4
Funnel plot
We used funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry.
Figure corresponds to Fig. 7A in main text
# funnel plot
Funnel_E <- funnel(mod_Efm, xlab = "lnRR", ylab = "Standard Error")# Funnel_EFig.S5 Funnel plot showing the standard error and residuals (lnRR) from the full model.
Egger’s regression and time-lag bias
We fitted effective sample size calculated as:
\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]
and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year.
VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster = dat_Efm$Study_ID,
r = 0.5)
dat_Efm$c_Year_published <- as.vector(scale(dat_Efm$Year_published, scale = F))
dat_Efm$sqrt_inv_e_n <- with(dat_Efm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
PB_MR_E <- rma.mv(lnRR_Sa, VCV_Efm, mods = ~1 + sqrt_inv_e_n + Learning_vs_memory +
Year_published + Type_assay + Type_reinforcement + EE_social + EE_exercise +
Age_stress_exposure, random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML",
test = "t", data = dat_Efm)
# PB_MR_E
estimates_PB_MR_E <- estimates.CI(PB_MR_E)
estimates_PB_MR_E## estimate mean lower upper
## 1 intrcpt -9.48981044 -50.89863818 31.91901730
## 2 sqrt_inv_e_n -0.06281042 -0.86857733 0.74295650
## 3 Learning_vs_memoryMemory -0.06858190 -0.17753636 0.04037255
## 4 Year_published 0.00467568 -0.01591231 0.02526367
## 5 Type_assayHabituation -0.55467640 -0.99644706 -0.11290575
## 6 Type_assayRecognition -0.11960483 -0.45315288 0.21394322
## 7 Type_reinforcementAversive 0.15536913 -0.15356140 0.46429965
## 8 Type_reinforcementNot applicable 0.32741009 -0.11084714 0.76566731
## 9 EE_socialSocial 0.04883160 -0.16642221 0.26408540
## 10 EE_exerciseNo exercise -0.03983198 -0.28097488 0.20131093
## 11 Age_stress_exposureAdult -0.16033127 -0.54043113 0.21976858
## 12 Age_stress_exposureEarly postnatal -0.05978991 -0.45195630 0.33237649
## 13 Age_stress_exposurePrenatal -0.13988907 -0.53943086 0.25965272
Leave-one-out analysis
We individually removed studies to determine if any singular studies have a disproportionate effect on the meta-analytic mean and CI.
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
d <- dat %>%
filter(Study_ID != levels(dat$Study_ID)[i])
VCV_Eb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Ea, V = VCV_Eb, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
levels(dat$Study_ID)[i], ])
}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0) {
df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
return(df)
}
# using dplyr to form data frame
MA_CVR_E <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
bind_rows %>%
mutate(left_out = levels(dat$Study_ID))
# saveRDS(MA_CVR_E,file = here('Rdata', 'MA_CVR_E.rds'))Figure corresponds to Fig. S3 in Supplementary material S2.
# MA_CVR_E <- readRDS(file = here('Rdata', 'MA_CVR_E.rds'))
# telling ggplot to stop reordering factors
MA_CVR_E$left_out <- as.factor(MA_CVR_E$left_out)
MA_CVR_E$left_out <- factor(MA_CVR_E$left_out, levels = MA_CVR_E$left_out)
# plotting
leaveoneout_E <- ggplot(MA_CVR_E) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_E0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_E0$ci.ub,
lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
leaveoneout_Edat$Study_ID <- as.integer(dat$Study_ID)Fig. S6 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.
Using SMD
We re-ran the main effect of environmental enrichment meta-analytic model using standardised mean difference (SMD) instead of lnRR.
mod_E0a <- rma.mv(yi = SMD_Ea, V = VCV_Ea, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_E0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -111.6284 223.2568 231.2568 241.3003 231.7220
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2512 0.5012 30 no Study_ID
## sigma^2.2 0.4247 0.6517 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 673.4722, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.7241 0.1295 5.5895 91 <.0001 0.4668 0.9814 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 8.673320e-01 3.223202e-01 5.450118e-01 1.749563e-09
Risk of Bias
Testing a model including if the study was randomised as a moderator and conducting a pair-wise contrast analysis.
mod_randomised_E <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~ROB_randomisation - 1,
random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_randomised_E)
r2_ml(mod_randomised_E)contra_mod_randomised_E <- contrast_fun(data = dat, response = lnRR_Ea, moderator = ROB_randomisation,
VCV = VCV_E)
res_table_mod_randomised_E <- get_estimate(model = mod_randomised_E, contra = contra_mod_randomised_E,
moderator = ROB_randomisation)
res_table_mod_randomised_E## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear 0.140 0.0488 0.232 0.00304 -0.276 0.556
## 2 Yes 0.227 0.139 0.315 0.00000158 -0.188 0.642
## 3 Unclear-Yes 0.0868 -0.0352 0.209 0.161 -0.337 0.511
Testing a model including if the study used blinding as a moderator and conducting a pair-wise contrast analysis.
mod_blinding_E <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~ROB_blinding - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_blinding_E)
r2_ml(mod_blinding_E)contra_mod_blinding_E <- contrast_fun(data = dat, response = lnRR_Ea, moderator = ROB_blinding,
VCV = VCV_E)
res_table_mod_blinding_E <- get_estimate(model = mod_blinding_E, contra = contra_mod_blinding_E,
moderator = ROB_blinding)
res_table_mod_blinding_E## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear 0.198 0.122 0.274 0.00000131 -0.226 0.622
## 2 Yes 0.148 0.00811 0.288 0.0384 -0.292 0.588
## 3 Unclear-Yes -0.0501 -0.208 0.108 0.530 -0.496 0.396
Main effect: stress
Meta-analysis
To quantify differences in individual performance in cognitive assays with stress, we used the logarithm of response ratio (lnRR) calculated as:
\[ \ln{\text{RR}_\text{stress}} = \ln\left( \frac{ M_\text{ES} + M_\text{CS}} {2} \right) - \ln\left( \frac {M_\text{EC} + M_\text{CC}} {2} \right) = \ln \left( {\frac {M_\text{ES} + M_\text{CS}} {M_\text{EC} + M_\text{CC}}} \right), \]
Variance was calculated as:
\[ \text{var}(\ln{\text{RR}_\text{stress}}) = \left( \frac{1} { M_\text{ES} + M_\text{CS} } \right)^2 \left( \frac{SD_\text{ES}^2}{N_\text{ES}} + \frac{SD_\text{CS}^2}{N_\text{CS}} \right) + \left( \frac {1} { M_\text{EC} + M_\text{CC} } \right)^2 \left(\frac{SD_\text{EC}^2}{N_\text{EC}} + \frac{SD_\text{CC}^2}{N_\text{CC}} \right) \] Variance was converted in to variance-covariance (VCV) matrix (see above in ‘Data organisation’) to control for non-independence in sampling variance from the same studies.
mod_S0 <- rma.mv(yi = lnRR_Sa, V = VCV_S, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_S0)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.1156 28.2312 36.2312 46.2747 36.6964
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0099 0.0993 30 no Study_ID
## sigma^2.2 0.0377 0.1941 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 946.9234, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1052 0.0337 -3.1217 91 0.0024 -0.1721 -0.0383 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.376124e-01 1.946168e-01 7.429955e-01 1.975455e-10
Figure corresponds to Fig. 3A in main text
orchard_plot(mod_S0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Fig. S7 Orchard plot showing meta-analytic mean and 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision.
Meta-regression: uni-moderator
To explain some of the unexplained variation in the main effect of stress model, we conducted a series of uni-moderator analyses. We calculated marginal R^2 for each moderator as well as conducted a series of pair-wise contrasts between moderator categories.
Any moderator categories with k < 5 were removed.
Learning vs Memory
Was the measured response interpreted as learning or memory?
mod_S1 <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Learning_vs_memory - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_S1)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.5281 29.0562 39.0562 51.5553 39.7705
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0094 0.0970 30 no Study_ID
## sigma^2.2 0.0388 0.1969 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 946.3930, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 5.0145, p-val = 0.0086
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning -0.1211 0.0476 -2.5423 90 0.0127 -0.2157
## Learning_vs_memoryMemory -0.0974 0.0380 -2.5648 90 0.0120 -0.1728
## ci.ub
## Learning_vs_memoryLearning -0.0265 *
## Learning_vs_memoryMemory -0.0219 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S1)## R2_marginal R2_coditional
## 0.002776034 1.000000000
Figure corresponds to Fig. 5A in main text
LvsM_S <- orchard_plot(mod_S1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_S Fig. S8a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_S1 <- contrast_fun(data = dat, response = lnRR_Sa, moderator = Learning_vs_memory,
VCV = VCV_S)
res_table_mod_S1 <- get_estimate(model = mod_S1, contra = contra_mod_S1, moderator = Learning_vs_memory)
res_table_mod_S1## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Learning -0.121 -0.216 -0.0265 0.0127 -0.567 0.325
## 2 Memory -0.0974 -0.173 -0.0219 0.0120 -0.540 0.345
## 3 Learning-Memory 0.0237 -0.0780 0.125 0.644 -0.424 0.472
Type of assay
The broad category of the type of assay used to measure learning or memory.
dat$Type_assay <- as.factor(dat$Type_assay)
VCV_S1 <- impute_covariance_matrix(vi = dat1$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S2 <- rma.mv(yi = lnRR_Sa, V = VCV_S1, mod = ~Type_assay - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat1)
summary(mod_S2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.8028 19.6055 31.6055 46.5374 32.6299
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0161 0.1271 30 no Study_ID
## sigma^2.2 0.0279 0.1671 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 723.4973, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 6.7053, p-val = 0.0004
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning -0.0981 0.0375 -2.6192 89 0.0104 -0.1725 -0.0237
## Type_assayHabituation -0.4615 0.1126 -4.0969 89 <.0001 -0.6853 -0.2377
## Type_assayRecognition -0.0534 0.0645 -0.8287 89 0.4095 -0.1816 0.0747
##
## Type_assayConditioning *
## Type_assayHabituation ***
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S2)## R2_marginal R2_coditional
## 0.1853359 1.0000000
Figure corresponds to Fig. 5B in main text
Learning_S <- orchard_plot(mod_S2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Habituation", "Recognition")),
annotations = "***") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_SFig. S8b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_S2 <- contrast_fun(data = dat1, response = lnRR_Sa, moderator = Type_assay,
VCV = VCV_S1)
res_table_mod_S2 <- get_estimate(model = mod_S2, contra = contra_mod_S2, moderator = Type_assay)
res_table_mod_S2## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Conditioning -0.0981 -0.173 -0.0237 1.04e-2 -0.522 0.326
## 2 Habituation -0.461 -0.685 -0.238 9.21e-5 -0.935 0.0119
## 3 Recognition -0.0534 -0.182 0.0747 4.10e-1 -0.490 0.383
## 4 Conditioning-Habituation -0.363 -0.584 -0.142 1.55e-3 -0.835 0.109
## 5 Conditioning-Recognition 0.0447 -0.0882 0.177 5.06e-1 -0.393 0.482
## 6 Habituation-Recognition 0.408 0.161 0.655 1.48e-3 -0.0768 0.893
Type of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
VCV_S2 <- impute_covariance_matrix(vi = dat2$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S3 <- rma.mv(yi = lnRR_Sa, V = VCV_S2, mod = ~Type_reinforcement - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat2)
summary(mod_S3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.8810 27.7621 39.7621 54.6939 40.7864
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0103 0.1016 30 no Study_ID
## sigma^2.2 0.0387 0.1966 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 920.8439, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.7942, p-val = 0.0130
##
## Model Results:
##
## estimate se tval df pval
## Type_reinforcementAppetitive -0.1846 0.0749 -2.4649 89 0.0156
## Type_reinforcementAversive -0.0730 0.0427 -1.7081 89 0.0911
## Type_reinforcementNot applicable -0.1172 0.0590 -1.9851 89 0.0502
## ci.lb ci.ub
## Type_reinforcementAppetitive -0.3334 -0.0358 *
## Type_reinforcementAversive -0.1579 0.0119 .
## Type_reinforcementNot applicable -0.2345 0.0001 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S3)## R2_marginal R2_coditional
## 0.0366428 1.0000000
Figure corresponds to Fig. 5C in main text
Reinforcement_S <- orchard_plot(mod_S3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_SFig. S8c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_S3 <- contrast_fun(data = dat2, response = lnRR_Sa, moderator = Type_reinforcement,
VCV = VCV_S2)
res_table_mod_S3 <- get_estimate(model = mod_S3, contra = contra_mod_S3, moderator = Type_reinforcement)
res_table_mod_S3## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Appetitive -0.185 -0.333 -3.58e-2 0.0156 -0.649 0.280
## 2 Aversive -0.0730 -0.158 1.19e-2 0.0911 -0.521 0.375
## 3 Not applicable -0.117 -0.235 1.12e-4 0.0502 -0.572 0.338
## 4 Appetitive-Aversive 0.112 -0.0592 2.82e-1 0.198 -0.360 0.583
## 5 Appetitive-Not applicable 0.0674 -0.119 2.54e-1 0.474 -0.410 0.545
## 6 Aversive-Not applicable -0.0442 -0.173 8.46e-2 0.497 -0.502 0.414
Age at stress
The age when individuals were exposed to stress.
mod_S4 <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Age_stress_exposure - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_S4)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -12.4083 24.8166 38.8166 56.1579 40.2166
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0048 0.0694 30 no Study_ID
## sigma^2.2 0.0392 0.1979 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 881.1229, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 4.3703, p-val = 0.0029
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent 0.0074 0.1159 0.0641 88 0.9490
## Age_stress_exposureAdult -0.2279 0.0622 -3.6664 88 0.0004
## Age_stress_exposureEarly postnatal -0.0561 0.0435 -1.2891 88 0.2007
## Age_stress_exposurePrenatal -0.1145 0.0743 -1.5404 88 0.1271
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.2228 0.2377
## Age_stress_exposureAdult -0.3514 -0.1044 ***
## Age_stress_exposureEarly postnatal -0.1425 0.0304
## Age_stress_exposurePrenatal -0.2621 0.0332
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S4)## R2_marginal R2_coditional
## 0.09987307 1.00000000
Figure corresponds to Fig. 5D in main text
Age_S <- orchard_plot(mod_S4, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_signif(comparisons = list(c("Adult", "Early postnatal")),
annotations = "*") +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_S Fig. S8d Orchard plot showing the group-wise means of the categorical variable ‘Age_stress_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_S4 <- contrast_fun(data = dat, response = lnRR_Sa, moderator = Age_stress_exposure,
VCV = VCV_S)
res_table_mod_S4 <- get_estimate(model = mod_S4, contra = contra_mod_S4, moderator = Age_stress_exposure)
res_table_mod_S4## # A tibble: 10 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adolescent 0.00743 -0.223 0.238 9.49e-1 -0.469 0.484
## 2 Adult -0.228 -0.351 -0.104 4.20e-4 -0.663 0.207
## 3 Early postnatal -0.0561 -0.142 0.0304 2.01e-1 -0.482 0.370
## 4 Prenatal -0.114 -0.262 0.0332 1.27e-1 -0.557 0.328
## 5 Adolescent-Adult -0.235 -0.497 0.0260 7.70e-2 -0.727 0.257
## 6 Adolescent-Early postna… -0.0635 -0.309 0.182 6.09e-1 -0.547 0.421
## 7 Adolescent-Prenatal -0.122 -0.395 0.152 3.78e-1 -0.620 0.377
## 8 Adult-Early postnatal 0.172 0.0211 0.323 2.60e-2 -0.271 0.615
## 9 Adult-Prenatal 0.113 -0.0791 0.306 2.45e-1 -0.346 0.573
## 10 Early postnatal-Prenatal -0.0584 -0.229 0.113 4.99e-1 -0.509 0.392
Type of stress
The type of stress used.
dat5 <- filter(dat, Type_stress_exposure %in% c("Restraint", "Noise", "Maternal separation",
"Combination"))
VCV_S3 <- impute_covariance_matrix(vi = dat5$lnRRV_S, cluster = dat5$Study_ID, r = 0.5)
mod_S5 <- rma.mv(yi = lnRR_Sa, V = VCV_S3, mod = ~Type_stress_exposure - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat5)
summary(mod_S5)##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -11.4138 22.8276 36.8276 53.5887 38.3618
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0115 0.1071 25 no Study_ID
## sigma^2.2 0.0401 0.2002 85 no ES_ID
## sigma^2.3 0.0000 0.0000 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 897.4553, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 2.8767, p-val = 0.0278
##
## Model Results:
##
## estimate se tval df pval
## Type_stress_exposureCombination -0.0500 0.0892 -0.5605 81 0.5767
## Type_stress_exposureMaternal separation -0.0539 0.0560 -0.9630 81 0.3384
## Type_stress_exposureNoise -0.1203 0.1036 -1.1608 81 0.2491
## Type_stress_exposureRestraint -0.2101 0.0704 -2.9863 81 0.0037
## ci.lb ci.ub
## Type_stress_exposureCombination -0.2274 0.1274
## Type_stress_exposureMaternal separation -0.1653 0.0575
## Type_stress_exposureNoise -0.3265 0.0859
## Type_stress_exposureRestraint -0.3501 -0.0701 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S5)## R2_marginal R2_coditional
## 0.07029554 1.00000000
Figure corresponds to Fig. 5E in main text
Stressor<- orchard_plot(mod_S5, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
StressorFig. S8e Orchard plot showing the group-wise means of the categorical variable ‘Type_stress_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_S5 <- contrast_fun(data = dat5, response = lnRR_Sa, moderator = Type_stress_exposure,
VCV = VCV_S3)
res_table_mod_S5 <- get_estimate(model = mod_S5, contra = contra_mod_S5, moderator = Type_stress_exposure)
res_table_mod_S5## # A tibble: 10 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Combination -0.0500 -0.227 0.127 0.577 -0.535 0.435
## 2 Maternal separation -0.0539 -0.165 0.0575 0.338 -0.519 0.411
## 3 Noise -0.120 -0.326 0.0859 0.249 -0.617 0.376
## 4 Restraint -0.210 -0.350 -0.0701 0.00373 -0.683 0.263
## 5 Combination-Maternal se… -0.00393 -0.213 0.206 0.970 -0.502 0.494
## 6 Combination-Noise -0.0703 -0.342 0.202 0.608 -0.598 0.457
## 7 Combination-Restraint -0.160 -0.386 0.0658 0.162 -0.665 0.345
## 8 Maternal separation-Noi… -0.0664 -0.301 0.168 0.575 -0.575 0.443
## 9 Maternal separation-Res… -0.156 -0.335 0.0227 0.0861 -0.642 0.330
## 10 Noise-Restraint -0.0898 -0.339 0.159 0.475 -0.606 0.426
Stess duration
Was the stress acute or chronic?
dat6 <- filter(dat, Stress_duration %in% c("Chronic", "Acute"))
VCV_S4 <- impute_covariance_matrix(vi = dat6$lnRRV_S, cluster = dat6$Study_ID, r = 0.5)
mod_S6 <- rma.mv(yi = lnRR_Sa, V = VCV_S4, mod = ~Stress_duration - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat6)
summary(mod_S6)##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.7898 27.5796 37.5796 49.9092 38.3204
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0062 0.0786 29 no Study_ID
## sigma^2.2 0.0391 0.1978 89 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 915.9393, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 6.6661, p-val = 0.0020
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute 0.0135 0.0659 0.2045 87 0.8384 -0.1176 0.1445
## Stress_durationChronic -0.1360 0.0373 -3.6456 87 0.0005 -0.2101 -0.0618
##
## Stress_durationAcute
## Stress_durationChronic ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S6)## R2_marginal R2_coditional
## 0.08245896 1.00000000
Figure corresponds to Fig. 5F in main text
Duration_S <- orchard_plot(mod_S6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Duration_S Fig. S8f Orchard plot showing the group-wise means of the categorical variable ‘Stress_duration’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_S6 <- contrast_fun(data = dat6, response = lnRR_Sa, moderator = Stress_duration,
VCV = VCV_S4)
res_table_mod_S6 <- get_estimate(model = mod_S6, contra = contra_mod_S6, moderator = Stress_duration)
res_table_mod_S6## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Acute 0.0135 -0.118 0.145 0.838 -0.429 0.456
## 2 Chronic -0.136 -0.210 -0.0618 0.000453 -0.565 0.294
## 3 Acute-Chronic -0.149 -0.300 0.00111 0.0517 -0.599 0.300
Publication bias & sensitivity analysis
Multi-moderator model
We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.
# selecting moderator levels with k >=5
dat_Sfm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
c("Appetitive", "Aversive", "Not applicable"), Type_stress_exposure %in%
c("Restraint", "Noise", "Maternal separation", "Combination"), Stress_duration %in%
c("Chronic", "Acute"))
# VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster =
# dat_Sfm$Study_ID, r = 0.5)
mod_Sfm <- rma.mv(yi = lnRR_Sa, V = VCV_Sfm, mod = ~Type_assay - 1 + Learning_vs_memory +
Type_reinforcement + Type_stress_exposure + Age_stress_exposure + Stress_duration,
random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Sfm)
# summary(mod_Sfm) r2_ml(mod_Sfm)
res_Sfm <- dredge(mod_Sfm, trace = 2)
saveRDS(res_Sfm, file = here("Rdata", "res_Sfm.rds"))
# also saving the full model and data
saveRDS(mod_Sfm, file = here("Rdata", "mod_Sfm.rds"))
saveRDS(dat_Sfm, file = here("Rdata", "dat_Sfm.rds"))The akaike weights for the top set of models with AIC < 6
dat_Sfm <- readRDS(file = here("Rdata", "dat_Sfm.rds"))
mod_Sfm <- readRDS(file = here("Rdata", "mod_Sfm.rds"))
res_Sfm <- readRDS(file = here("Rdata", "res_Sfm.rds"))
res_Sfm2 <- subset(res_Sfm, delta <= 6, recalc.weights = FALSE)
importance(res_Sfm2)## Type_assay Stress_duration Type_stress_exposure
## Sum of weights: 0.91 0.88 0.22
## N containing models: 7 6 2
## Learning_vs_memory Age_stress_exposure Type_reinforcement
## Sum of weights: 0.16 0.04 0.04
## N containing models: 2 1 1
Funnel plot
We used funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry.
Figure corresponds to Fig. 7B in main text.
# funnel plot
Funnel_S <- funnel(mod_Sfm, xlab = "lnRR", ylab = "Standard Error")# Funnel_SFig.S9 Funnel plot showing the standard error and residuals (lnRR) from the full model.
Egger’s regression and time-lag bias
We fitted effective sample size calculated as:
\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]
and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year.
# calculating inv effective sample size for use in full meta-regression
VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster = dat_Sfm$Study_ID,
r = 0.5)
dat_Sfm$sqrt_inv_e_n <- with(dat_Sfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
# time lag bias and eggers regression
dat_Sfm$c_Year_published <- as.vector(scale(dat_Sfm$Year_published, scale = F))
PB_MR_S <- rma.mv(lnRR_Sa, VCV_Sfm, mods = ~1 + sqrt_inv_e_n + c_Year_published +
Type_assay + Learning_vs_memory + Type_reinforcement + Type_stress_exposure +
Age_stress_exposure + Stress_duration, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), method = "REML", test = "t", data = dat_Sfm, control = list(optimizer = "optim",
optmethod = "Nelder-Mead"))
# PB_MR_S
estimates_PB_MR_S <- estimates.CI(PB_MR_S)
estimates_PB_MR_S## estimate mean lower upper
## 1 intrcpt 0.208757092 -0.47118467 0.88869885
## 2 sqrt_inv_e_n 0.109877719 -0.90636737 1.12612281
## 3 c_Year_published -0.006425821 -0.02804904 0.01519740
## 4 Type_assayHabituation -0.663428076 -1.09349872 -0.23335743
## 5 Type_assayRecognition -0.069789046 -0.41156145 0.27198335
## 6 Learning_vs_memoryMemory -0.076712961 -0.19562432 0.04219839
## 7 Type_reinforcementAversive 0.039563987 -0.13761334 0.21674132
## 8 Type_reinforcementNot applicable 0.189481532 -0.16955283 0.54851589
## 9 Type_stress_exposureMaternal separation 0.233030567 -0.40478886 0.87084999
## 10 Type_stress_exposureNoise -0.020700600 -0.76667312 0.72527192
## 11 Type_stress_exposureRestraint -0.067490308 -0.43366600 0.29868539
## 12 Age_stress_exposureAdult -0.129875123 -0.51040393 0.25065368
## 13 Age_stress_exposureEarly postnatal -0.193461086 -0.89857848 0.51165631
## 14 Age_stress_exposurePrenatal -0.146406300 -0.57377820 0.28096560
## 15 Stress_durationChronic -0.395711964 -0.65566557 -0.13575836
Leave-one-out sensitivity analysis
We individually removed studies to determine if any singular studies have a disproportionate effect on the meta-analytic mean and CI.
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
d <- dat %>%
filter(Study_ID != levels(dat$Study_ID)[i])
VCV_Sb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Sa, V = VCV_Sb, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
levels(dat$Study_ID)[i], ])
}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_S0) {
df <- data.frame(est = mod_S0$b, lower = mod_S0$ci.lb, upper = mod_S0$ci.ub)
return(df)
}
# using dplyr to form data frame
MA_CVR_S <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
bind_rows %>%
mutate(left_out = levels(dat$Study_ID))
# saveRDS(MA_CVR_S,file = here('Rdata', 'MA_CVR_S.rds'))Figure corresponds to Fig. S4 in Supplementary material S2
# MA_CVR_S <- readRDS(file = here('Rdata', 'MA_CVR_S.rds'))
# telling ggplot to stop reordering factors
MA_CVR_S$left_out <- as.factor(MA_CVR_S$left_out)
MA_CVR_S$left_out <- factor(MA_CVR_S$left_out, levels = MA_CVR_S$left_out)
# plotting
leaveoneout_S <- ggplot(MA_CVR_S) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_S0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_S0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_S0$ci.ub,
lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
leaveoneout_Sdat$Study_ID <- as.integer(dat$Study_ID)Fig. S10 Leave-one-study-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.
Using SMD
We re-ran the main effect of stress meta-analytic model using standardised mean difference (SMD) instead of lnRR.
mod_S0a <- rma.mv(yi = SMD_Sa, V = VCV_Sa, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_S0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -120.2817 240.5634 248.5634 258.6068 249.0285
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3874 0.6224 30 no Study_ID
## sigma^2.2 0.4969 0.7049 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 818.5592, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.4256 0.1498 -2.8403 91 0.0056 -0.7232 -0.1279 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 8.959540e-01 3.924790e-01 5.034750e-01 1.438454e-09
Risk of Bias
Including if the study was randomised as a moderator and conducting a pair-wise contrast.
mod_randomised_S <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~ROB_randomisation - 1,
random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_randomised_S)
r2_ml(mod_randomised_S)contra_mod_randomised_S <- contrast_fun(data = dat, response = lnRR_Sa, moderator = ROB_randomisation,
VCV = VCV_S)
res_table_mod_randomised_S <- get_estimate(model = mod_randomised_S, contra = contra_mod_randomised_S,
moderator = ROB_randomisation)
res_table_mod_randomised_S## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear -0.0672 -0.164 0.0291 0.169 -0.509 0.374
## 2 Yes -0.139 -0.230 -0.0480 0.00317 -0.580 0.301
## 3 Unclear-Yes -0.0719 -0.204 0.0607 0.284 -0.523 0.379
Including if the study used blinding and conducting a pair-wise contrast.
mod_blinding_S <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~ROB_blinding - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_blinding_S)
r2_ml(mod_blinding_S)contra_mod_blinding_S <- contrast_fun(data = dat, response = lnRR_Sa, moderator = ROB_blinding,
VCV = VCV_S)
res_table_mod_blinding_S <- get_estimate(model = mod_blinding_S, contra = contra_mod_blinding_S,
moderator = ROB_blinding)
res_table_mod_blinding_S## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear -0.119 -0.196 -0.0419 0.00288 -0.564 0.326
## 2 Yes -0.0567 -0.202 0.0890 0.442 -0.518 0.405
## 3 Unclear-Yes 0.0625 -0.102 0.227 0.454 -0.406 0.531
Interaction: enrichment x stress
Meta-analysis
To quantify differences in individual performance in cognitive assay with the interaction of environmental enrichment and stress, we used the logarithm of response ratio (lnRR) calculated as:
\[ \ln{\text{RR}_\text{interaction}} = \left(\ln M_\text{ES} - \ln M_\text{CS} \right) - \left( \ln M_\text{EC} - \ln M_\text{CC} \right) = \ln \left( \frac {M_\text{ES} / M_\text{CS} } {M_\text{EC} / \ M_\text{CC}} \right) \]
Variance was calculated as:
\[ \text{var}(\ln{\text{RR}_\text{interaction}}) = \frac {SD_\text{ES}^2} { N_\text{ES} M_\text{ES}^2 } + \frac{SD_\text{EC}^2}{N_\text{EC}M_\text{EC}^2} + \frac {SD_\text{CS}^2}{N_\text{CS} M_\text{CS}^2} + \frac{SD_\text{CC}^2}{N_\text{CC} M_\text{CC}^2}. \]
Variance was converted into variance-covariance (VCV) matrix (see above in ‘Data organisation’) to control for non-independence in sampling variance from the same studies.
mod_ES0 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_ES0)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.8178 81.6355 89.6355 99.6790 90.1006
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0316 0.1777 30 no Study_ID
## sigma^2.2 0.0229 0.1513 92 no ES_ID
## sigma^2.3 0.0030 0.0544 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 303.2179, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1229 0.0596 2.0605 91 0.0422 0.0044 0.2414 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 0.81703063 0.44913873 0.32576306 0.04212884
Figure corresponds to Fig. 3A in main text
orchard_plot(mod_ES0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Fig. S11 Orchard plot showing meta-analytic mean and 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Meta-regression: uni-moderator
To explain some of the unexplained variation in the interaction of environmental enrichment and stress, we conducted a series of uni-moderator analyses. We calculated marginal R^2 for each moderator as well as conducted a series of pair-wise contrasts between moderator categories.
Any moderator categories with k < 5 were removed.
Learning vs Memory
Was the response interpreted as learning or memory?
mod_ES1 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Learning_vs_memory - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_ES1)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.4769 80.9539 90.9539 103.4529 91.6682
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0292 0.1708 30 no Study_ID
## sigma^2.2 0.0232 0.1524 92 no ES_ID
## sigma^2.3 0.0034 0.0582 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 299.1854, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 2.9219, p-val = 0.0590
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.1744 0.0722 2.4166 90 0.0177 0.0310
## Learning_vs_memoryMemory 0.1057 0.0619 1.7065 90 0.0914 -0.0174
## ci.ub
## Learning_vs_memoryLearning 0.3178 *
## Learning_vs_memoryMemory 0.2288 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES1)## R2_marginal R2_coditional
## 0.0197648 0.9405080
Figure corresponds to Fig. 6A in main text
LvsM_ES <- orchard_plot(mod_ES1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 15), # change font sizes
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 10))
LvsM_ES # axis.text.x = element_text(size = 10), # change font sizes
# axis.text.y = element_text(size = 10),
# legend.title = element_text(size = 7),
# legend.text = element_text(size = 7)) Fig. S12a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_ES1 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Learning_vs_memory,
VCV = VCV_ES)
res_table_mod_ES1 <- get_estimate(model = mod_ES1, contra = contra_mod_ES1, moderator = Learning_vs_memory)
res_table_mod_ES1## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Learning 0.174 0.0310 0.318 0.0177 -0.316 0.665
## 2 Memory 0.106 -0.0174 0.229 0.0914 -0.379 0.591
## 3 Learning-Memory -0.0687 -0.176 0.0382 0.205 -0.550 0.413
Type of assay
The broad category of the type of assay used to measure learning or memory.
VCV_ES1 <- impute_covariance_matrix(vi = dat1$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
mod_ES2 <- rma.mv(yi = lnRR_ESa, V = VCV_ES1, mod = ~Type_assay - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat1)
summary(mod_ES2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0874 78.1748 90.1748 105.1066 91.1992
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0370 0.1924 30 no Study_ID
## sigma^2.2 0.0192 0.1386 92 no ES_ID
## sigma^2.3 0.0018 0.0422 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.9385, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1062, p-val = 0.0305
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning 0.1525 0.0589 2.5877 89 0.0113 0.0354 0.2696
## Type_assayHabituation 0.1990 0.1415 1.4070 89 0.1629 -0.0820 0.4801
## Type_assayRecognition -0.0048 0.0800 -0.0606 89 0.9518 -0.1637 0.1541
##
## Type_assayConditioning *
## Type_assayHabituation
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES2)## R2_marginal R2_coditional
## 0.05775809 0.97105550
Figure corresponds to Fig. 6B in main text
Learning_ES <- orchard_plot(mod_ES2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Conditioning", "Recognition")),
annotations = "*") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 15), # change font sizes
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 10))
Learning_ESFig. S12b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_ES2 <- contrast_fun(data = dat1, response = lnRR_ESa, moderator = Type_assay,
VCV = VCV_ES1)
res_table_mod_ES2 <- get_estimate(model = mod_ES2, contra = contra_mod_ES2, moderator = Type_assay)
res_table_mod_ES2## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Conditioning 0.153 0.0354 0.270 0.0113 -0.340 0.645
## 2 Habituation 0.199 -0.0820 0.480 0.163 -0.356 0.754
## 3 Recognition -0.00485 -0.164 0.154 0.952 -0.509 0.499
## 4 Conditioning-Habituation 0.0465 -0.219 0.312 0.728 -0.501 0.594
## 5 Conditioning-Recognition -0.157 -0.304 -0.0105 0.0360 -0.658 0.343
## 6 Habituation-Recognition -0.204 -0.488 0.0801 0.157 -0.760 0.353
Type of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
VCV_ES2 <- impute_covariance_matrix(vi = dat2$lnRRV_ES, cluster = dat2$Study_ID,
r = 0.5)
mod_ES3 <- rma.mv(yi = lnRR_ESa, V = VCV_ES2, mod = ~Type_reinforcement - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat2)
summary(mod_ES3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0604 78.1208 90.1208 105.0526 91.1452
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0382 0.1954 30 no Study_ID
## sigma^2.2 0.0189 0.1377 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.4724, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.2547, p-val = 0.0254
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.1007 0.1075 0.9366 89 0.3515 -0.1129
## Type_reinforcementAversive 0.1573 0.0569 2.7618 89 0.0070 0.0441
## Type_reinforcementNot applicable 0.0147 0.0702 0.2101 89 0.8341 -0.1247
## ci.ub
## Type_reinforcementAppetitive 0.3143
## Type_reinforcementAversive 0.2704 **
## Type_reinforcementNot applicable 0.1541
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES3)## R2_marginal R2_coditional
## 0.0586952 1.0000000
Figure corresponds to Fig. 6C in main text
Reinforcement_ES <- orchard_plot(mod_ES3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
geom_signif(comparisons = list(c("Aversive", "Not applicable")),
annotations = "*") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 15), # change font sizes
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 10))
Reinforcement_ES Fig. S12c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_ES3 <- contrast_fun(data = dat2, response = lnRR_ESa, moderator = Type_reinforcement,
VCV = VCV_ES2)
res_table_mod_ES3 <- get_estimate(model = mod_ES3, contra = contra_mod_ES3, moderator = Type_reinforcement)
res_table_mod_ES3## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Appetitive 0.101 -0.113 0.314 0.351 -0.420 0.621
## 2 Aversive 0.157 0.0441 0.270 0.00698 -0.331 0.645
## 3 Not applicable 0.0147 -0.125 0.154 0.834 -0.480 0.510
## 4 Appetitive-Aversive 0.0566 -0.182 0.295 0.639 -0.475 0.588
## 5 Appetitive-Not applicable -0.0859 -0.333 0.161 0.492 -0.621 0.450
## 6 Aversive-Not applicable -0.143 -0.278 -0.00744 0.0389 -0.636 0.351
Age at enrichment
The age when individuals were exposed to environmental enrichment.
VCV_ES3 <- impute_covariance_matrix(vi = dat3$lnRRV_ES, cluster = dat3$Study_ID,
r = 0.5)
mod_ES4 <- rma.mv(yi = lnRR_ESa, V = VCV_ES3, mod = ~Age_EE_exposure - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat3)
summary(mod_ES4)##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -36.6407 73.2813 83.2813 95.5531 84.0313
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0332 0.1822 29 no Study_ID
## sigma^2.2 0.0207 0.1439 88 no ES_ID
## sigma^2.3 0.0007 0.0265 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 288.6493, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 3.1308, p-val = 0.0487
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Age_EE_exposureAdolescent 0.1284 0.0583 2.2006 86 0.0304 0.0124
## Age_EE_exposureAdult 0.1247 0.0921 1.3538 86 0.1793 -0.0584
## ci.ub
## Age_EE_exposureAdolescent 0.2444 *
## Age_EE_exposureAdult 0.3077
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES4)## R2_marginal R2_coditional
## 0.0000400928 0.9871095822
Figure corresponds to Fig. 6D in main text
Age_enrichment_ES <- orchard_plot(mod_ES4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 15), # change font sizes
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 10))
Age_enrichment_ESFig. S12d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_ES4 <- contrast_fun(data = dat3, response = lnRR_ESa, moderator = Age_EE_exposure,
VCV = VCV_ES3)
res_table_mod_ES4 <- get_estimate(model = mod_ES4, contra = contra_mod_ES4, moderator = Age_EE_exposure)
res_table_mod_ES4## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adolescent 0.128 0.0124 0.244 0.0304 -0.350 0.607
## 2 Adult 0.125 -0.0584 0.308 0.179 -0.375 0.624
## 3 Adolescent-Adult -0.00373 -0.213 0.205 0.972 -0.513 0.506
Age at stress
The age when individuals were exposed to stress.
mod_ES7 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Age_stress_exposure - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_ES7)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0917 78.1834 92.1834 109.5247 93.5834
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0295 0.1718 30 no Study_ID
## sigma^2.2 0.0232 0.1522 92 no ES_ID
## sigma^2.3 0.0091 0.0954 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 286.4225, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 1.5932, p-val = 0.1832
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent -0.0137 0.1762 -0.0775 88 0.9384
## Age_stress_exposureAdult 0.1677 0.1140 1.4708 88 0.1449
## Age_stress_exposureEarly postnatal 0.1067 0.0920 1.1602 88 0.2491
## Age_stress_exposurePrenatal 0.3179 0.1311 2.4254 88 0.0173
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.3639 0.3366
## Age_stress_exposureAdult -0.0589 0.3942
## Age_stress_exposureEarly postnatal -0.0761 0.2896
## Age_stress_exposurePrenatal 0.0574 0.5784 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES7)## R2_marginal R2_coditional
## 0.09276574 0.86630830
Figure corresponds to Fig. 6E in main text
Age_stress_ES<-orchard_plot(mod_ES7, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 15), # change font sizes
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 10))
Age_stress_ESFig. S12e Orchard plot showing the group-wise means of the categorical variable ‘Age_stress_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_ES7 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Age_stress_exposure,
VCV = VCV_ES)
res_table_mod_ES7 <- get_estimate(model = mod_ES7, contra = contra_mod_ES7, moderator = Age_stress_exposure)
res_table_mod_ES7## # A tibble: 10 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adolescent -0.0137 -0.364 0.337 0.938 -0.619 0.592
## 2 Adult 0.168 -0.0589 0.394 0.145 -0.376 0.711
## 3 Early postnatal 0.107 -0.0761 0.290 0.249 -0.420 0.634
## 4 Prenatal 0.318 0.0574 0.578 0.0173 -0.241 0.876
## 5 Adolescent-Adult 0.181 -0.236 0.598 0.390 -0.465 0.828
## 6 Adolescent-Early postna… 0.120 -0.275 0.515 0.546 -0.512 0.753
## 7 Adolescent-Prenatal 0.332 -0.105 0.768 0.135 -0.328 0.991
## 8 Adult-Early postnatal -0.0609 -0.280 0.158 0.581 -0.601 0.479
## 9 Adult-Prenatal 0.150 -0.132 0.432 0.292 -0.419 0.719
## 10 Early postnatal-Prenatal 0.211 -0.0462 0.469 0.107 -0.346 0.768
Order to treatment exposure
The order in which individuals were exposed to enrichment and stress.
mod_ES10 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Exposure_order - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_ES10)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0408 78.0817 90.0817 105.0135 91.1061
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0316 0.1777 30 no Study_ID
## sigma^2.2 0.0227 0.1507 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 292.2561, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1546, p-val = 0.0287
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Exposure_orderConcurrently 0.1492 0.1208 1.2351 89 0.2201 -0.0909
## Exposure_orderEnrichment first -0.1782 0.1659 -1.0744 89 0.2856 -0.5079
## Exposure_orderStress first 0.1370 0.0526 2.6046 89 0.0108 0.0325
## ci.ub
## Exposure_orderConcurrently 0.3893
## Exposure_orderEnrichment first 0.1514
## Exposure_orderStress first 0.2414 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES10)## R2_marginal R2_coditional
## 0.1027207 1.0000000
Figure corresponds to Fig. 6F in main text
Order_ES <- orchard_plot(mod_ES10, mod = "Exposure_order", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 15), # change font sizes
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 10))
Order_ES Fig. S12f Orchard plot showing the group-wise means of the categorical variable ‘Exposure_order’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_ES10 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Exposure_order,
VCV = VCV_ES)
res_table_mod_ES10 <- get_estimate(model = mod_ES10, contra = contra_mod_ES10, moderator = Exposure_order)
res_table_mod_ES10## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Concurrently 0.149 -0.0909 0.389 0.220 -0.372 0.671
## 2 Enrichment first -0.178 -0.508 0.151 0.286 -0.747 0.390
## 3 Stress first 0.137 0.0325 0.241 0.0108 -0.338 0.612
## 4 Concurrently-Enrichment … -0.327 -0.735 0.0803 0.114 -0.944 0.289
## 5 Concurrently-Stress first -0.0123 -0.274 0.250 0.926 -0.544 0.520
## 6 Enrichment first-Stress … 0.315 -0.0306 0.661 0.0735 -0.263 0.893
Exercise enrichment
Did enrichment involve the addition of apparatus for voluntary exercise in the cage (e.g., a wheel or treadmill)?
mod_ES5 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~EE_exercise - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_ES5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.4968 80.9935 90.9935 103.4926 91.7078
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0325 0.1803 30 no Study_ID
## sigma^2.2 0.0230 0.1517 92 no ES_ID
## sigma^2.3 0.0065 0.0805 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 297.3270, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 1.8401, p-val = 0.1647
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1051 0.0780 1.3474 90 0.1812 -0.0499 0.2602
## EE_exerciseNo exercise 0.1687 0.0967 1.7448 90 0.0844 -0.0234 0.3609
##
## EE_exerciseExercise
## EE_exerciseNo exercise .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES5)## R2_marginal R2_coditional
## 0.01474459 0.89712469
Figure corresponds to Fig. 6G in main text
Exercise_ES <- orchard_plot(mod_ES5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 15), # change font sizes
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 10))
Exercise_ESFig. S12g Orchard plot showing the group-wise means of the categorical variable ‘EE_exercise’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_ES5 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = EE_exercise,
VCV = VCV_ES)
res_table_mod_ES5 <- get_estimate(model = mod_ES5, contra = contra_mod_ES5, moderator = EE_exercise)
res_table_mod_ES5## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Exercise 0.105 -0.0499 0.260 0.181 -0.413 0.624
## 2 No exercise 0.169 -0.0234 0.361 0.0844 -0.362 0.699
## 3 Exercise-No exercise 0.0636 -0.138 0.265 0.532 -0.470 0.598
Type of stress
The type of stressor used.
VCV_ES5 <- impute_covariance_matrix(vi = dat5$lnRRV_ES, cluster = dat5$Study_ID,
r = 0.5)
mod_ES8 <- rma.mv(yi = lnRR_ESa, V = VCV_ES5, mod = ~Type_stress_exposure - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat5)
summary(mod_ES8)##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -34.4046 68.8091 82.8091 99.5703 84.3434
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0426 0.2064 25 no Study_ID
## sigma^2.2 0.0232 0.1524 85 no ES_ID
## sigma^2.3 0.0104 0.1021 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 281.9708, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 0.5137, p-val = 0.7258
##
## Model Results:
##
## estimate se tval df pval
## Type_stress_exposureCombination 0.1111 0.1458 0.7618 81 0.4484
## Type_stress_exposureMaternal separation 0.1185 0.1099 1.0784 81 0.2840
## Type_stress_exposureNoise 0.1651 0.1795 0.9198 81 0.3604
## Type_stress_exposureRestraint 0.1374 0.1252 1.0978 81 0.2755
## ci.lb ci.ub
## Type_stress_exposureCombination -0.1790 0.4011
## Type_stress_exposureMaternal separation -0.1001 0.3370
## Type_stress_exposureNoise -0.1920 0.5221
## Type_stress_exposureRestraint -0.1116 0.3865
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES8)## R2_marginal R2_coditional
## 0.004455703 0.863910284
Figure corresponds to Fig. 6I in main text
Stressor_ES <- orchard_plot(mod_ES8, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 15), # change font sizes
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 10))
Stressor_ES Fig. S12I Orchard plot showing the group-wise means of the categorical variable ‘Type_stress_exposure’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_ES8 <- contrast_fun(data = dat5, response = lnRR_ESa, moderator = Type_stress_exposure,
VCV = VCV_ES5)
res_table_mod_ES8 <- get_estimate(model = mod_ES8, contra = contra_mod_ES8, moderator = Type_stress_exposure)
res_table_mod_ES8## # A tibble: 10 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Combination 0.111 -0.179 0.401 0.448 -0.510 0.732
## 2 Maternal separation 0.118 -0.100 0.337 0.284 -0.473 0.710
## 3 Noise 0.165 -0.192 0.522 0.360 -0.490 0.820
## 4 Restraint 0.137 -0.112 0.386 0.276 -0.466 0.741
## 5 Combination-Maternal se… 0.00741 -0.299 0.313 0.962 -0.621 0.636
## 6 Combination-Noise 0.0540 -0.352 0.460 0.792 -0.629 0.737
## 7 Combination-Restraint 0.0263 -0.300 0.353 0.873 -0.613 0.665
## 8 Maternal separation-Noi… 0.0466 -0.312 0.405 0.797 -0.609 0.703
## 9 Maternal separation-Res… 0.0189 -0.246 0.284 0.887 -0.591 0.629
## 10 Noise-Restraint -0.0276 -0.403 0.348 0.884 -0.693 0.638
Stress duration
Was the stress acute or chronic?
VCV_ES6 <- impute_covariance_matrix(vi = dat6$lnRRV_ES, cluster = dat6$Study_ID,
r = 0.5)
mod_ES9 <- rma.mv(yi = lnRR_ESa, V = VCV_ES6, mod = ~Stress_duration - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat6)
summary(mod_ES9)##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -35.9010 71.8020 81.8020 94.1315 82.5427
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0133 0.1155 29 no Study_ID
## sigma^2.2 0.0260 0.1611 89 no ES_ID
## sigma^2.3 0.0080 0.0895 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 278.4877, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 4.3877, p-val = 0.0153
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute -0.0367 0.0930 -0.3946 87 0.6941 -0.2216 0.1482
## Stress_durationChronic 0.1854 0.0732 2.5347 87 0.0130 0.0400 0.3308
##
## Stress_durationAcute
## Stress_durationChronic *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES9)## R2_marginal R2_coditional
## 0.1598311 0.8575918
Figure corresponds to Fig. 6J in main text
Duration_ES<- orchard_plot(mod_ES9, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_signif(comparisons = list(c("Acute", "Chronic")),
annotations = "*") +
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 15), # change font sizes
axis.text.y = element_text(size = 15),
legend.title = element_text(size = 10))
Duration_ESFig. S12j Orchard plot showing the group-wise means of the categorical variable ‘Stress_duration’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts
contra_mod_ES9 <- contrast_fun(data = dat6, response = lnRR_ESa, moderator = Stress_duration,
VCV = VCV_ES6)
res_table_mod_ES9 <- get_estimate(model = mod_ES9, contra = contra_mod_ES9, moderator = Stress_duration)
res_table_mod_ES9## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Acute -0.0367 -0.222 0.148 0.694 -0.507 0.433
## 2 Chronic 0.185 0.0400 0.331 0.0130 -0.271 0.642
## 3 Acute-Chronic 0.222 0.0381 0.406 0.0186 -0.248 0.692
Publication bias & sensitivity analysis
Multi-moderator model
We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.
dat_ESfm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
c("Appetitive", "Aversive", "Not applicable"), EE_social %in% c("Social",
"Non-social"), Age_EE_exposure %in% c("Adult", "Adolescent"), Type_stress_exposure %in%
c("Restraint", "Noise", "Maternal separation", "Combination"), Stress_duration %in%
c("Chronic", "Acute"), Age_stress_exposure %in% c("Prenatal", "Early postnatal",
"Adult"))
# VCV_ESfm <- impute_covariance_matrix(vi = dat_ESfm$lnRRV_ES, cluster =
# dat_ESfm$Study_ID, r = 0.5)
mod_ESfm <- rma.mv(yi = lnRR_Sa, V = VCV_ESfm, mod = ~Type_assay - 1 + Learning_vs_memory +
Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure + Type_stress_exposure +
Age_stress_exposure + Stress_duration + Exposure_order, random = list(~1 | Study_ID,
~1 | ES_ID, ~1 | Strain), test = "t", data = dat_ESfm)
# summary(mod_ESfm) r2_ml(mod_ESfm)
res_ESfm <- dredge(mod_ESfm, trace = 2)
saveRDS(res_ESfm, file = here("Rdata", "res_ESfm.rds"))
# also saving the full model and data
saveRDS(mod_ESfm, file = here("Rdata", "mod_ESfm.rds"))
saveRDS(dat_ESfm, file = here("Rdata", "dat_ESfm.rds"))Akaike weights for the top set of models with AIC < 6
dat_ESfm <- readRDS(file = here("Rdata", "dat_ESfm.rds"))
mod_ESfm <- readRDS(file = here("Rdata", "mod_ESfm.rds"))
res_ESfm <- readRDS(file = here("Rdata", "res_ESfm.rds"))
res_ESfm2 <- subset(res_ESfm, delta <= 6, recalc.weights = FALSE)
importance(res_ESfm2)## Type_assay Stress_duration Age_EE_exposure EE_exercise
## Sum of weights: 0.65 0.63 0.26 0.11
## N containing models: 19 18 7 6
## Learning_vs_memory EE_social Age_stress_exposure
## Sum of weights: 0.08 0.08 0.05
## N containing models: 4 4 2
## Type_stress_exposure Exposure_order
## Sum of weights: 0.03 0.02
## N containing models: 2 1
Funnel plot
We used funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry.
Figure corresponds to Fig. 7C in main text
Funnel_ES <- funnel(mod_ESfm, xlab = "lnRR", ylab = "Standard Error")# Funnel_ESFig.S13 Funnel plot showing the standard error and residuals (lnRR) from the full model.
Egger’s regression and time-lag bias
We fitted effective sample size calculated as:
\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]
and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year.
# calculating inv effective sample size for use in full meta-regression
dat_ESfm$sqrt_inv_e_n <- with(dat_ESfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
dat_ESfm$c_Year_published <- as.vector(scale(dat_ESfm$Year_published, scale = F))
VCV_ESfm <- impute_covariance_matrix(vi = dat_ESfm$lnRRV_ES, cluster = dat_ESfm$Study_ID,
r = 0.5)
# time lag bias and Egger's regression.
PB_MR_ES <- rma.mv(lnRR_ESa, VCV_ESfm, mods = ~1 + sqrt_inv_e_n + c_Year_published +
Type_assay + Learning_vs_memory + Type_reinforcement + Type_stress_exposure +
Age_stress_exposure + EE_social + EE_exercise + Stress_duration + Exposure_order,
random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", test = "t",
data = dat_ESfm, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
estimates_PB_MR_ES <- estimates.CI(PB_MR_ES)
estimates_PB_MR_ES## estimate mean lower upper
## 1 intrcpt 0.104890890 -1.32422406 1.53400584
## 2 sqrt_inv_e_n -0.380617946 -2.34738544 1.58614955
## 3 c_Year_published -0.003614793 -0.03254341 0.02531382
## 4 Type_assayHabituation 0.225637592 -0.28343491 0.73471009
## 5 Type_assayRecognition 0.014042786 -0.37390986 0.40199543
## 6 Learning_vs_memoryMemory 0.027907733 -0.10104526 0.15686073
## 7 Type_reinforcementAversive 0.074150779 -0.20556264 0.35386420
## 8 Type_reinforcementNot applicable -0.180164458 -0.64972276 0.28939384
## 9 Type_stress_exposureMS -0.670185502 -1.72077062 0.38039961
## 10 Type_stress_exposureNoise -0.307675466 -1.41676992 0.80141899
## 11 Type_stress_exposureRestraint -0.159808179 -0.87246493 0.55284857
## 12 Age_stress_exposureEarly postnatal 0.472151143 -0.39118506 1.33548735
## 13 Age_stress_exposurePrenatal 0.199706727 -0.39069154 0.79010499
## 14 EE_socialSocial 0.134840904 -0.25767695 0.52735875
## 15 EE_exerciseNo exercise 0.127498329 -0.18382738 0.43882404
## 16 Stress_durationChronic 0.449057668 0.05381018 0.84430516
## 17 Exposure_orderEnrichment first -0.132255288 -0.74038190 0.47587133
## 18 Exposure_orderStress first -0.057436396 -0.54737789 0.43250509
Leave-one-out analysis
We individually removed studies to determine if any singular studies have a disproportionate effect on the meta-analytic mean and CI.
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
d <- dat %>%
filter(Study_ID != levels(dat$Study_ID)[i])
VCV_ESb <- impute_covariance_matrix(vi = d$lnRRV_ES, cluster = d$Study_ID, r = 0.5)
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_ESa, V = VCV_ESb, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
levels(dat$Study_ID)[i], ])
}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_ES0) {
df <- data.frame(est = mod_ES0$b, lower = mod_ES0$ci.lb, upper = mod_ES0$ci.ub)
return(df)
}
# using dplyr to form data frame
MA_CVR_ES <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
bind_rows %>%
mutate(left_out = levels(dat$Study_ID))
# saveRDS(MA_CVR_ES, ,file = here('Rdata', 'MA_CVR_ES.rds'))Figure corresponds to Fig. S5 in Supplementary material S2
# MA_CVR_ES<- readRDS(here('Rdata', 'MA_CVR_ES.rds'))
# telling ggplot to stop reordering factors
MA_CVR_ES$left_out <- as.factor(MA_CVR_ES$left_out)
MA_CVR_ES$left_out <- factor(MA_CVR_ES$left_out, levels = MA_CVR_ES$left_out)
# plotting
leaveoneout_ES <- ggplot(MA_CVR_ES) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_ES0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_ES0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_ES0$ci.ub,
lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
leaveoneout_ESdat$Study_ID <- as.integer(dat$Study_ID)Fig. S14 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.
Using SMD
We re-ran the interaction of environmental enrichment x stress meta-analytic model using standardised mean difference (SMD) instead of lnRR.
mod_ES0a <- rma.mv(yi = SMD_ESa, V = VCV_ESa, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_ES0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -126.0571 252.1141 260.1141 270.1576 260.5793
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4653 0.6821 30 no Study_ID
## sigma^2.2 0.3698 0.6081 92 no ES_ID
## sigma^2.3 0.0000 0.0002 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 257.4673, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.6880 0.1763 3.9017 91 0.0002 0.3377 1.0382 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 6.657130e-01 3.709364e-01 2.947766e-01 3.351907e-08
Risk of Bias
Including if the study was randomised as a moderator and conducting a pair-wise contrast.
mod_randomised_ES <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~ROB_randomisation -
1, random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_randomised_ES)
r2_ml(mod_randomised_ES)contra_mod_randomised_ES <- contrast_fun(data = dat, response = lnRR_ESa, moderator = ROB_randomisation,
VCV = VCV_ES)
res_table_mod_randomised_ES <- get_estimate(model = mod_randomised_ES, contra = contra_mod_randomised_ES,
moderator = ROB_randomisation)
res_table_mod_randomised_ES## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear 0.0467 -0.106 0.200 0.546 -0.434 0.527
## 2 Yes 0.200 0.0503 0.350 0.00946 -0.279 0.680
## 3 Unclear-Yes 0.154 -0.0200 0.327 0.0821 -0.333 0.641
Including if the study used blinding and conducting a pair-wise contrast
mod_blinding_ES <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~ROB_blinding - 1, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)
summary(mod_blinding_ES)
r2_ml(mod_blinding_ES)contra_mod_blinding_ES <- contrast_fun(data = dat, response = lnRR_ESa, moderator = ROB_blinding,
VCV = VCV_ES)
res_table_mod_blinding_ES <- get_estimate(model = mod_blinding_ES, contra = contra_mod_blinding_ES,
moderator = ROB_blinding)
res_table_mod_blinding_ES## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Unclear 0.158 -0.000309 0.316 0.0504 -0.364 0.680
## 2 Yes 0.0444 -0.178 0.267 0.693 -0.500 0.589
## 3 Unclear-Yes -0.114 -0.343 0.115 0.327 -0.661 0.434
‘Pairwise’ effect sizes
We conducted five ‘meta-analyses on ’pair-wise’ comparisons using lnRR between the treatments and control:
\[
\ln{\text{RR}} = \ln (M_{\text{EC}}/M_{\text{CC}}),
\ln (M_{\text{CS}}/M_{\text{CC}}),
\ln (M_{\text{ES}}/M_{\text{CC}}),
\ln (M_{\text{ES}}/M_{\text{CS}})
\ln (M_{\text{ES}}/M_{\text{EC}}),
\]
Enrichment relative to control
\[ \ln{\text{RR}} = \ln (M_{\text{EC}}/M_{\text{CC}}) \]
VCV_E20 <- impute_covariance_matrix(vi = dat$lnRRV_E2, cluster = dat$Study_ID, r = 0.5)
mod_E20 <- rma.mv(yi = lnRR_E2a, V = VCV_E20, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
summary(mod_E20)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.3235 14.6470 22.6470 32.6904 23.1121
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0037 0.0611 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0281 0.1675 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 475.8327, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1066 0.0291 3.6655 91 0.0004 0.0489 0.1644 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E20)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.610789e-01 1.011724e-01 2.607945e-09 7.599065e-01
Figure corresponds to Fig. 3B in main text
orchard_plot(mod_E20, mod = "Int", xlab = "lnRR")Fig. S15a Orchard plot showing meta-analytic mean and 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Stress relative to control
\[ \ln{\text{RR}} = \ln (M_{\text{CS}}/M_{\text{CC}}) \]
VCV_S20 <- impute_covariance_matrix(vi = dat$lnRRV_S2, cluster = dat$Study_ID, r = 0.5)
mod_S20 <- rma.mv(yi = lnRR_S2a, V = VCV_S20, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat)
summary(mod_S20)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -52.3561 104.7122 112.7122 122.7557 113.1773
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0323 0.1797 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0798 0.2824 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 1003.0694, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1846 0.0522 -3.5360 91 0.0006 -0.2882 -0.0809 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S20)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 9.489604e-01 2.734220e-01 9.879730e-10 6.755384e-01
Figure corresponds to Fig. 3B in main text
orchard_plot(mod_S20, mod = "Int", xlab = "lnRR")Fig. S15b Orchard plot showing meta-analytic mean and 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Enrichment + stress relative to control
\[ \ln{\text{RR}} = \ln (M_{\text{ES}}/M_{\text{CC}}) \]
VCV_ES20 <- impute_covariance_matrix(vi = dat$lnRRV_ES2, cluster = dat$Study_ID,
r = 0.5)
mod_ES20 <- rma.mv(yi = lnRR_ES2a, V = VCV_ES20, random = list(~1 | Study_ID, ~1 |
Strain, ~1 | ES_ID), test = "t", data = dat)
summary(mod_ES20)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.3625 20.7250 28.7250 38.7684 29.1901
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0039 0.0625 30 no Study_ID
## sigma^2.2 0.0014 0.0377 6 no Strain
## sigma^2.3 0.0227 0.1508 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 402.2656, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.0801 0.0389 2.0594 91 0.0423 0.0028 0.1573 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES20)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 0.81513970 0.11355557 0.04132363 0.66026050
Figure corresponds to Fig. 3B in main text
orchard_plot(mod_ES20, mod = "Int", xlab = "lnRR")Fig. S15c Orchard plot showing meta-analytic mean and 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Enrichment + stress relative to stress
\[ \ln{\text{RR}} = \ln (M_{\text{ES}}/M_{\text{CS}}) \]
VCV_E30 <- impute_covariance_matrix(vi = dat$lnRRV_E3, cluster = dat$Study_ID, r = 0.5)
mod_E30 <- rma.mv(yi = lnRR_E3a, V = VCV_E30, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat)
summary(mod_E30)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -46.3447 92.6895 100.6895 110.7329 101.1546
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0235 0.1532 30 no Study_ID
## sigma^2.2 0.0263 0.1623 6 no Strain
## sigma^2.3 0.0543 0.2331 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 790.0249, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.2465 0.1011 2.4389 91 0.0167 0.0457 0.4472 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E30)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 0.9456920 0.2132063 0.2391809 0.4933048
Figure corresponds to Fig. 3B in main text
orchard_plot(mod_E30, mod = "Int", xlab = "lnRR")Fig. S15d Orchard plot showing meta-analytic mean and 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Enrichment + stress relative to enrichment
\[ \ln{\text{RR}} = \ln (M_{\text{ES}}/M_{\text{EC}}) \]
VCV_S30 <- impute_covariance_matrix(vi = dat$lnRRV_S3, cluster = dat$Study_ID, r = 0.5)
mod_S30 <- rma.mv(yi = lnRR_S3a, V = VCV_S30, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
summary(mod_S30)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -6.8788 13.7576 21.7576 31.8011 22.2228
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0000 0.0000 30 no Study_ID
## sigma^2.2 0.0009 0.0292 6 no Strain
## sigma^2.3 0.0276 0.1662 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 540.3522, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.0087 0.0342 -0.2552 91 0.7992 -0.0766 0.0592
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S30)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.415428e-01 5.192043e-09 2.515933e-02 8.163835e-01
Figure correspinds to Fig. 3B in main text
orchard_plot(mod_S30, mod = "Int", xlab = "lnRR")Fig. S15e Orchard plot showing meta-analytic mean and 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Figures
Panel of ‘focal’ ES and ‘pairwise’ ES orchard plots
Figure corresponds to Fig.3 in main text
mod_list1 <- list(mod_E0, mod_S0, mod_ES0)
mod_res1 <- lapply(mod_list1, function(x) mod_results(x, mod = "Int"))
merged1 <- submerge(mod_res1[[3]], mod_res1[[2]], mod_res1[[1]], mix = T)
merged1$mod_table$name <- factor(merged1$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
merged1$data$moderator <- factor(merged1$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
orchard1 <- orchard_plot(merged1, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 4, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling +
scale_colour_manual(values = c("#00AEEF","#00A651","#ED1C24"))+ # change colours
scale_fill_manual(values=c("#00AEEF","#00A651","#ED1C24"))+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10)) mod_list2 <- list(mod_S30, mod_E30, mod_ES20, mod_S20, mod_E20) #rearranged the order so that it matches intext results
mod_res2 <- lapply(mod_list2, function(x) mod_results(x, mod = "Int"))
merged2 <- submerge(mod_res2[[1]], mod_res2[[2]], mod_res2[[3]], mod_res2[[4]], mod_res2[[5]], mix = T)
merged2$mod_table$name <- factor(merged2$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
merged2$data$moderator <- factor(merged2$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
orchard2 <- orchard_plot(merged2, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 4, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
scale_colour_manual(values = c("#7B81BE","#D7DF23","#F37158","#75CBF2","#97D2B4"))+ # change colours
scale_fill_manual(values=c("#7B81BE","#D7DF23","#F37158","#75CBF2","#97D2B4"))+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10)) p1 <- orchard1 + orchard2 + plot_annotation(tag_levels = "A")
p1# saved as PDF: 6 x 15 inches.Fig. S16 Orchard plots showing the eight meta-analytic models. (A) main effects of enrichment and stress, and interaction of environmental enrichment and stress, (B) ‘pairwise’ comparisons between treatments and controls. Thick black lines = 95% CI, thin black lines = 95 % PI. Individual effect sizes are shown as circles scaled proportionately to their precision. k is number of effect sizes.
Panel of meta-regressions
Environmental enrichment
Figure corresponds to Fig. 4 in main text
# Enrichment
E_mod <- (LvsM_E + Learning_E + Reinforcement_E)/(Age_E + Exercise_E + Social_E) +
plot_annotation(tag_levels = "A")
E_mod# saved as pdf 10 x 15 inchesFig. S17 Orchard plots showing the six different meta-regressions on the main effect of environmental enrichment on learning and memory. (A) learning versus memory response, (B) the type of assay, (C) the type of reinforcement, (D) the age at environmental enrichment, (E) type of manipulation of exercise during enrichment, (F) manipulation of the social environment during enrichment. Individual efect sizes are shown as circles scaled proportionately to their precision. k is number of effect sizes.
Stress
Figure corresponds to Fig. 5 in main text
S_mod <- (LvsM_S + Learning_S + Reinforcement_S)/(Age_S + Stressor + Duration_S) +
plot_annotation(tag_levels = "A")
S_mod# saved as pdf 10 x 15 inchesFig. S18 Orchard plots showing the six different meta-regressions on the main effect of stress on learning and memory. (A) learning versus memory response, (B) the type of assay, (C) the type of reinforcement, (D) the age at stress, (E) the type of stressor, (F) chronic or acute stress. Individual effect sizes are shown as circles scaled proportionately to their precision. k is number of effect sizes.
Interaction
Figure corresponds to Fig. 6 in main text
ES_mod <- plot_grid(LvsM_ES, Learning_ES, Reinforcement_ES, Age_enrichment_ES, Age_stress_ES,
Order_ES, Exercise_ES, Social_ES, Stressor_ES, Duration_ES, labels = "AUTO",
ncol = 5)
ES_mod# saved as 10 x 25 inchesFig. S19 Orchard plots showing the 10 different meta-regressions of moderators on the interaction between environmental enrichment and stress learning and memory. (A) learning versus memory response, (B) the type of assay, (C) the type of reinforcement used, (D) the age at environmental enrichment, (E) the age at stress, (F) the order of treatment exposure, (G) if enrichment involved a manipulation of exercise, (H) manipulation of the social environment during enrichment, (I) the type of stressor, (F) stress was chronic or acute. Individual effect sizes are shown as circles scaled proportionately to their precision. k is number of effect sizes.
Panel of funnel plots
Figure corresponds to Fig. 7 in main text
# EE
pdf(NULL)
dev.control(displaylist = "enable")
par(mar = c(4, 4, 0.1, 0))
A <- funnel(mod_Sfm, xlab = "Residuals (lnRR)", ylab = "Standard Error", xlim = c(-2,
2), ylim = c(0, 1.05))
A <- recordPlot()
invisible(dev.off())
# Stress
pdf(NULL)
dev.control(displaylist = "enable")
par(mar = c(4, 4, 0.1, 0))
B <- funnel(mod_Sfm, xlab = "Residuals (lnRR)", ylab = "Standard Error", xlim = c(-2,
2), ylim = c(0, 1.05))
B <- recordPlot()
invisible(dev.off())
# Interaction
pdf(NULL)
dev.control(displaylist = "enable")
par(mar = c(4, 4, 0.1, 0))
C <- funnel(mod_ESfm, xlab = "Residuals (lnRR)", ylab = "Standard Error", xlim = c(-2,
2), ylim = c(0, 1.05))
C <- recordPlot()
invisible(dev.off())
# putting together
ggdraw(A) + ggdraw(B) + ggdraw(C) + plot_annotation(tag_levels = "A")
# png(file = here('figs', 'Fig7_Funnels.png'))
# dev.off()knitr::include_graphics(here("figs", "funnels.png"))Fig. S20 Funnel plots of the standard error and residuals (lnRR) from the full models. (A) environmental enrichment main effect, (B) stress main effect, (C) environmental enrichment x stress interaction.
Software and package versions
sessionInfo() %>%
pander()R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8
attached base packages: grid, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: formatR(v.1.11), pander(v.0.6.4), gridGraphics(v.0.5-1), png(v.0.1-7), cowplot(v.1.1.1), ggthemr(v.1.1.0), ggalluvial(v.0.12.3), visdat(v.0.5.3), ggsignif(v.0.6.3), networkD3(v.0.4), GoodmanKruskal(v.0.0.3), patchwork(v.1.1.1), MuMIn(v.1.43.17), orchaRd(v.0.0.0.9000), clubSandwich(v.0.5.3), metafor(v.3.0-2), Matrix(v.1.3-3), here(v.1.0.1), forcats(v.0.5.1), stringr(v.1.4.0), dplyr(v.1.0.6), purrr(v.0.3.4), readr(v.1.4.0), tidyr(v.1.1.3), tibble(v.3.1.2), ggplot2(v.3.3.5) and tidyverse(v.1.3.1)
loaded via a namespace (and not attached): nlme(v.3.1-152), fs(v.1.5.0), lubridate(v.1.7.10), RColorBrewer(v.1.1-2), httr(v.1.4.2), rprojroot(v.2.0.2), tools(v.4.1.0), backports(v.1.2.1), utf8(v.1.2.1), R6(v.2.5.0), vipor(v.0.4.5), DBI(v.1.1.1), colorspace(v.2.0-1), withr(v.2.4.2), tidyselect(v.1.1.1), compiler(v.4.1.0), cli(v.2.5.0), rvest(v.1.0.0), pacman(v.0.5.1), xml2(v.1.3.2), sandwich(v.3.0-1), labeling(v.0.4.2), bookdown(v.0.22), scales(v.1.1.1), digest(v.0.6.27), rmarkdown(v.2.8), pkgconfig(v.2.0.3), htmltools(v.0.5.1.1), highr(v.0.9), dbplyr(v.2.1.1), htmlwidgets(v.1.5.3), rlang(v.0.4.11), readxl(v.1.3.1), rstudioapi(v.0.13), farver(v.2.1.0), generics(v.0.1.0), zoo(v.1.8-9), jsonlite(v.1.7.2), magrittr(v.2.0.1), ggbeeswarm(v.0.6.0), Rcpp(v.1.0.6), munsell(v.0.5.0), fansi(v.0.4.2), lifecycle(v.1.0.0), stringi(v.1.6.2), yaml(v.2.2.1), mathjaxr(v.1.4-0), crayon(v.1.4.1), lattice(v.0.20-44), haven(v.2.4.1), hms(v.1.1.0), knitr(v.1.33), pillar(v.1.6.1), igraph(v.1.2.6), stats4(v.4.1.0), reprex(v.2.0.0), glue(v.1.4.2), evaluate(v.0.14), modelr(v.0.1.8), vctrs(v.0.3.8), rmdformats(v.1.0.2), cellranger(v.1.1.0), gtable(v.0.3.0), assertthat(v.0.2.1), xfun(v.0.23), broom(v.0.7.6), beeswarm(v.0.4.0) and ellipsis(v.0.3.2)
Social enrichment
Did enrichment involve adding more conspecifics in the cage? (Note that we did not include studies that only provided social enrichment but no abiotic enrichment.)
Figure corresponds to Fig. 4F in main text
Fig. S4f Orchard plot showing the group-wise means of the categorical variable ‘EE_social’ with their 95% confidence interval (thick black line) and 95% prediction interval (thin black line). Individual points show observed effect sizes scaled proportionately to their precision. k is number of effect sizes.
Contrasts